/* Include files */
#include "math.h"
#include "mex.h"

/* Include to sort : http://stackoverflow.com/questions/1787996/c-library-function-to-do-sort */
#include <stdio.h>
#include <stdlib.h>


/* Checks whether inputs are matrices */
#define IS_REAL_2D_FULL_DOUBLE(P) (!mxIsComplex(P) && \
mxGetNumberOfDimensions(P) == 2 && !mxIsSparse(P) && mxIsDouble(P))
#define IS_REAL_SCALAR(P) (IS_REAL_2D_FULL_DOUBLE(P) && mxGetNumberOfElements(P) == 1)

/*************************************************/
/* SORTING */

/* Sort elements with qsort */
/*
int comp (const void * elem1, const void * elem2) 
{
	int f = *((int*)elem1);
    int s = *((int*)elem2);
    if (f > s) return  1;
    if (f < s) return -1;
    return 0;
}
*/

/* Sort index by utility */
double *base_arr;
static int compPointer (const void *a, const void *b)
{
	int aa = *((int *) a), bb = *((int *) b);
	if (base_arr[aa - 1] < base_arr[bb - 1])
		return 1;
	if (base_arr[aa - 1] == base_arr[bb - 1])
		return 0;
	if (base_arr[aa - 1] > base_arr[bb - 1])
		return -1;
	return 0;
}

/*************************************************/
/* Check Aoffer is inputing by returning in Bset */
/*
void printMatEl(double *y, double *z, int m, int n)
{	
	int A = m; int a = 0;
	int B = n; int b = 0;
	for (a = 0; a < A; a++)
		for (b = 0; b < B; b++)
		{
			z[a + A*b] = y[a + A*b];
		}
}
*/
/*************************************************/

void assignVectorLen(double *y, double *z, int M, int N, int n)
{	
	int m = 0;
	for (m = 0; m < M; m++)
	{
		z[m] = y[m + M*n];
	}
}


/*************************************************/
/* MAIN */
/*************************************************/

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	/* Define output matrix */
	#define rankB	plhs[0]

	/* Define input matrices */
	#define Baccept 	prhs[0]

	int a, b;

	/* Check the number of inputs is correct (Aoffer,Baccept,QbarA) */
	if(nrhs < 1)
		mexErrMsgTxt("Too few input arguments.");
	else if(nrhs > 1)
		mexErrMsgTxt("Too many input arguments.");

	if(!IS_REAL_2D_FULL_DOUBLE(Baccept))
		mexErrMsgTxt("Baccept must be real 2D full double array.");
	
	/* Dimensions and pointers for inputs */
	int A = mxGetM(Baccept);
	int B = mxGetN(Baccept);
	int nb = A*B;

	double *B_pt = mxGetPr(Baccept);

	/* Output matrix */
	rankB = mxCreateDoubleMatrix(A,B,mxREAL);
	double *Rank_pt = mxGetPr(rankB);

	/* Uninitialized containers */
	mxArray *SortUtil = mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetData(SortUtil, mxMalloc(sizeof(double)*1*A));
	double *Util_pt = mxGetPr(SortUtil);

	/* Generate index to use for sorting */
	int	*idx = malloc (sizeof (int) * A);
	for (a = 0; a < A; a++)
	{
		idx[a] = a + 1;
/*		printf("index %d\n", idx[b]); */
	}

	/* Assignment test */
/*	printMatEl(A_pt,Bset_pt,A,B); */

/***************************************************************/
/* Begin Sorting */
/***************************************************************/

/* Assign pointer to vector of "offer" utilities */
	for (b = 0; b < B; b++)
	{
		
	/* Sort the index based on utility values (output = idx) */

		assignVectorLen(B_pt,Util_pt,A,B,b);
		base_arr = (double *)mxGetData(SortUtil);
		int *idxSort = idx;
		qsort (idxSort, A, sizeof(int), compPointer);


	/* Store the index result in the rankA matrix for each firm */
		for (a = 0; a < A; a++)
	    {
/*			Rank_pt[a + A*b] = idxSort[a]; */
	    	if (B_pt[(idxSort[a] - 1) + A*b] >= 0)
	    		Rank_pt[(idxSort[a] - 1) + A*b] = a + 1;
	    	else
	    		Rank_pt[(idxSort[a] - 1) + A*b] = 0;
	    }	

	}

	/* Check input */
/*  mexCallMATLAB(0,NULL,1, &Aoffer, "disp"); */

	mxDestroyArray(SortUtil);
	free(idx);

/*	return 0; */

}
